pacman::p_load(sf, spNetwork, tmap, tidyverse)Hands-on Ex 3
Network Constrained Spatial Point Pattern Analysis
1. Overview
Network-constrained spatial point pattern analysis (NetSPAA) focuses on studying events that occur along networks, such as roads and rivers, rather than across open spaces. Unlike traditional spatial analysis, which assumes events can happen anywhere and measures distances in straight lines, this method addresses phenomena restricted to networks, like traffic accidents or store locations. It is especially important for its accurate in detecting event clusters that occur along these networks.
The spNetwork package is a crucial geospatial R tool for conducting NetSPAA, specifically for:
- Network Kernel Density Estimation (NKDE), and
- Network-based G-function and K-function analysis.
2. Importing & Transforming Data
This study analyses the spatial distribution of childcare centers in the Punggol Planning Area. We will import and work with two geospatial datasets:
Punggol_St: A line feature dataset representing the road network within the Punggol Planning Area.Punggol_CC: A point feature dataset representing the locations of childcare centers within the Punggol Planning Area.
After importing the road network data, it is essential to review the dataframe, verify the coordinate reference system (CRS) is correct, and check for any duplicated geometries.
network <- st_read(dsn="data/geospatial",
layer="Punggol_St")Reading layer `Punggol_St' from data source
`/Users/stephentay/stephentay/ISSS626-Geospatial-Analytics/Hands-on_Ex/Hands-on_Ex03/data/geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 2642 features and 2 fields
Geometry type: LINESTRING
Dimension: XY
Bounding box: xmin: 34038.56 ymin: 40941.11 xmax: 38882.85 ymax: 44801.27
Projected CRS: SVY21 / Singapore TMst_crs(network)Coordinate Reference System:
User input: SVY21 / Singapore TM
wkt:
PROJCRS["SVY21 / Singapore TM",
BASEGEOGCRS["SVY21",
DATUM["SVY21",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4757]],
CONVERSION["Singapore Transverse Mercator",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",1.36666666666667,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",103.833333333333,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",1,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",28001.642,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",38744.572,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["northing (N)",north,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["easting (E)",east,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Cadastre, engineering survey, topographic mapping."],
AREA["Singapore - onshore and offshore."],
BBOX[1.13,103.59,1.47,104.07]],
ID["EPSG",3414]]networkSimple feature collection with 2642 features and 2 fields
Geometry type: LINESTRING
Dimension: XY
Bounding box: xmin: 34038.56 ymin: 40941.11 xmax: 38882.85 ymax: 44801.27
Projected CRS: SVY21 / Singapore TM
First 10 features:
LINK_ID ST_NAME geometry
1 116130894 PUNGGOL RD LINESTRING (36546.89 44574....
2 116130897 PONGGOL TWENTY-FOURTH AVE LINESTRING (36546.89 44574....
3 116130901 PONGGOL SEVENTEENTH AVE LINESTRING (36012.73 44154....
4 116130902 PONGGOL SEVENTEENTH AVE LINESTRING (36062.81 44197....
5 116130907 PUNGGOL CENTRAL LINESTRING (36131.85 42755....
6 116130908 PUNGGOL RD LINESTRING (36112.93 42752....
7 116130909 PUNGGOL CENTRAL LINESTRING (36127.4 42744.5...
8 116130910 PUNGGOL FLD LINESTRING (35994.98 42428....
9 116130911 PUNGGOL FLD LINESTRING (35984.97 42407....
10 116130912 EDGEFIELD PLNS LINESTRING (36200.87 42219....any(duplicated(network))[1] FALSEAfter importing the point data, it is essential to review the dataframe, verify the CRS is correct, and check for any duplicated geometries. Since the point geometries are in 3-dimensional coordinates, they must be transformed to 2-dimensional coordinates using the st_zm() method.
childcare <- st_read(dsn="data/geospatial",
layer="Punggol_CC")Reading layer `Punggol_CC' from data source
`/Users/stephentay/stephentay/ISSS626-Geospatial-Analytics/Hands-on_Ex/Hands-on_Ex03/data/geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 61 features and 1 field
Geometry type: POINT
Dimension: XYZ
Bounding box: xmin: 34423.98 ymin: 41503.6 xmax: 37619.47 ymax: 44685.77
z_range: zmin: 0 zmax: 0
Projected CRS: SVY21 / Singapore TMst_crs(childcare)Coordinate Reference System:
User input: SVY21 / Singapore TM
wkt:
PROJCRS["SVY21 / Singapore TM",
BASEGEOGCRS["SVY21",
DATUM["SVY21",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4757]],
CONVERSION["Singapore Transverse Mercator",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",1.36666666666667,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",103.833333333333,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",1,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",28001.642,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",38744.572,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["northing (N)",north,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["easting (E)",east,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Cadastre, engineering survey, topographic mapping."],
AREA["Singapore - onshore and offshore."],
BBOX[1.13,103.59,1.47,104.07]],
ID["EPSG",3414]]childcareSimple feature collection with 61 features and 1 field
Geometry type: POINT
Dimension: XYZ
Bounding box: xmin: 34423.98 ymin: 41503.6 xmax: 37619.47 ymax: 44685.77
z_range: zmin: 0 zmax: 0
Projected CRS: SVY21 / Singapore TM
First 10 features:
Name geometry
1 kml_10 POINT Z (36173.81 42550.33 0)
2 kml_99 POINT Z (36479.56 42405.21 0)
3 kml_100 POINT Z (36618.72 41989.13 0)
4 kml_101 POINT Z (36285.37 42261.42 0)
5 kml_122 POINT Z (35414.54 42625.1 0)
6 kml_161 POINT Z (36545.16 42580.09 0)
7 kml_172 POINT Z (35289.44 44083.57 0)
8 kml_188 POINT Z (36520.56 42844.74 0)
9 kml_205 POINT Z (36924.01 41503.6 0)
10 kml_222 POINT Z (37141.76 42326.36 0)any(duplicated(childcare))[1] FALSEchildcare$geometry <- st_zm(childcare$geometry)
childcareSimple feature collection with 61 features and 1 field
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 34423.98 ymin: 41503.6 xmax: 37619.47 ymax: 44685.77
Projected CRS: SVY21 / Singapore TM
First 10 features:
Name geometry
1 kml_10 POINT (36173.81 42550.33)
2 kml_99 POINT (36479.56 42405.21)
3 kml_100 POINT (36618.72 41989.13)
4 kml_101 POINT (36285.37 42261.42)
5 kml_122 POINT (35414.54 42625.1)
6 kml_161 POINT (36545.16 42580.09)
7 kml_172 POINT (35289.44 44083.57)
8 kml_188 POINT (36520.56 42844.74)
9 kml_205 POINT (36924.01 41503.6)
10 kml_222 POINT (37141.76 42326.36)3. GeoVisualisation of Childcare Centers
Before starting the analysis, it is good practice to visualize the geospatial data.
plot(st_geometry(network))
plot(childcare,add=T,col='red',pch = 19)
tmap_mode('view')tmap mode set to interactive viewingtm_shape(childcare) +
tm_dots() +
tm_shape(network) +
tm_lines()tmap_mode('plot')tmap mode set to plotting4. Network Kernel Density Estimation (NKDE) Analysis
4.1 Prepare lixels
Before computing NKDE, the SpatialLines object must be segmented into lixels with a specified minimum distance. The lixel length (lx_length) is set to 700m, and the minimum length (mindist) is set to 350m. After segmentation, if the final lixel is shorter than the minimum distance, it is merged with the previous lixel.
lixels <- lixelize_lines(network,
lx_length = 700,
mindist = 350)4.2 Compute lines’ center point
Next, lines_center() from the spNetwork package will be used to generate a SpatialPointsDataFrame with points positioned at the center of each line, based on the line’s length, as shown in the code below.
samples <- lines_center(lixels) 4.3 Perform NKDE
We are ready to computer the NKDE. There are 3 methods in computing NKDE:
The simple method, proposed by Xie and Yan (2008), differs from classical KDE in two key ways:
- Events are snapped to a network.
- Distances between sampling points and events are measured along the network (i.e., not using Euclidean distances).
Let’s examine the NKDE plot using the simple method below!
simple_densities <- nkde(lines = network, events = childcare,
w = rep(1, nrow(childcare)), samples = samples,
kernel_name = "quartic", bw = 300,
div= "bw", method = "simple",
digits = 1, tol = 1,
grid_shape = c(1,1), max_depth = 8,
agg = 5, sparse = TRUE,
verbose = FALSE)samples$simple_density <- simple_densities*1000
lixels$simple_density <- simple_densities*1000
tmap_mode('view')tmap mode set to interactive viewingtm_shape(lixels) +
tm_lines(col="simple_density") +
tm_shape(childcare) +
tm_dots()tmap_mode('plot')tmap mode set to plottingThe simple method has a limitation: the event’s mass is multiplied by the number of edges at an intersection, causing the kernel not to integrate to 1. To address this, Okabe and Sugihara (2012) proposed the discontinuous NKDE, which was further extended by Sugihara, Satoh, and Okabe (2010) for networks with cycles shorter than the bandwidth.
See the NKDE plot using the discontinuous method below!
discontinuous_densities <- nkde(lines = network, events = childcare,
w = rep(1, nrow(childcare)), samples = samples,
kernel_name = "quartic", bw = 300,
div= "bw", method = "discontinuous",
digits = 1, tol = 1,
grid_shape = c(1,1), max_depth = 8,
agg = 5, sparse = TRUE,
verbose = FALSE)samples$discontinuous_density <- discontinuous_densities*1000
lixels$discontinuous_density <- discontinuous_densities*1000
tmap_mode('view')tmap mode set to interactive viewingtm_shape(lixels) +
tm_lines(col="discontinuous_density") +
tm_shape(childcare) +
tm_dots()tmap_mode('plot')tmap mode set to plottingDiscontinuous NKDE can be counterintuitive, resulting in sharp differences between density values across the network, particularly problematic in networks with many closely spaced intersections. The continuous method adjusts NKDE values at intersections to ensure the kernel integrates to one, applying a backward correction to maintain continuity in density values. However, this method is recursive and computationally time-consuming.
See the NKDE plot using the continuous method below!
continuous_densities <- nkde(lines = network, events = childcare,
w = rep(1, nrow(childcare)), samples = samples,
kernel_name = "quartic", bw = 300,
div= "bw", method = "continuous",
digits = 1, tol = 1,
grid_shape = c(1,1), max_depth = 8,
agg = 5, sparse = TRUE,
verbose = FALSE)samples$continuous_density <- continuous_densities*1000
lixels$continuous_density <- continuous_densities*1000
tmap_mode('view')tmap mode set to interactive viewingtm_shape(lixels) +
tm_lines(col="continuous_density") +
tm_shape(childcare) +
tm_dots()tmap_mode('plot')tmap mode set to plotting5. Network-constrained G- and K-function analysis
A test for complete spatial randomness (CSR) can be performed using the network-constrained G- and K-functions.
The test hypotheses are as follows:
- H0: The distribution of childcare services in Punggol are randomly distributed along the road network.
- H1: The distribution of childcare services in Punggol are not randomly distributed along the road network.
Alpha is set at 0.05, with 1000 Monte-carlo simulations.
set.seed(2024)
kfun_childcare <- kfunctions(network,
childcare,
start = 0,
end = 1000,
step = 50,
width = 50,
nsim = 1000,
resolution = 50,
verbose = FALSE,
conf_int = 0.05)The blue line represents the empirical network K-function of the childcare centers in Punggol, while the gray area shows the results of 1,000 simulations within the 2.5% - 97.5% confidence interval. Since the blue line falls below the lower boundary of the envelope between the 200-400m distance range, it indicates a regular pattern of childcare center locations along the Punggol road network.
kfun_childcare$plotk
The blue line represents the empirical network G-function of the childcare centers in Punggol. Since the blue line falls below the lower boundary of the envelope in the 100-200m distance range, it suggests a regular pattern of childcare center locations along the Punggol road network.
kfun_childcare$plotg
6. Summary
Network-constrained spatial point pattern analysis revealed a regular distribution of childcare centers along the Punggol road network.